home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / quadpack / dqagie.f < prev    next >
Text File  |  1996-07-19  |  18KB  |  459 lines
  1.       SUBROUTINE DQAGIE(F,BOUND,INF,EPSABS,EPSREL,LIMIT,RESULT,ABSERR,
  2.      *   NEVAL,IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST)
  3. C***BEGIN PROLOGUE  DQAGIE
  4. C***DATE WRITTEN   800101   (YYMMDD)
  5. C***REVISION DATE  830518   (YYMMDD)
  6. C***CATEGORY NO.  H2A3A1,H2A4A1
  7. C***KEYWORDS  AUTOMATIC INTEGRATOR, INFINITE INTERVALS,
  8. C             GENERAL-PURPOSE, TRANSFORMATION, EXTRAPOLATION,
  9. C             GLOBALLY ADAPTIVE
  10. C***AUTHOR  PIESSENS,ROBERT,APPL. MATH & PROGR. DIV - K.U.LEUVEN
  11. C           DE DONCKER,ELISE,APPL. MATH & PROGR. DIV - K.U.LEUVEN
  12. C***PURPOSE  THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
  13. C            INTEGRAL   I = INTEGRAL OF F OVER (BOUND,+INFINITY)
  14. C            OR I = INTEGRAL OF F OVER (-INFINITY,BOUND)
  15. C            OR I = INTEGRAL OF F OVER (-INFINITY,+INFINITY),
  16. C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
  17. C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I))
  18. C***DESCRIPTION
  19. C
  20. C INTEGRATION OVER INFINITE INTERVALS
  21. C STANDARD FORTRAN SUBROUTINE
  22. C
  23. C            F      - DOUBLE PRECISION
  24. C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  25. C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  26. C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  27. C
  28. C            BOUND  - DOUBLE PRECISION
  29. C                     FINITE BOUND OF INTEGRATION RANGE
  30. C                     (HAS NO MEANING IF INTERVAL IS DOUBLY-INFINITE)
  31. C
  32. C            INF    - DOUBLE PRECISION
  33. C                     INDICATING THE KIND OF INTEGRATION RANGE INVOLVED
  34. C                     INF = 1 CORRESPONDS TO  (BOUND,+INFINITY),
  35. C                     INF = -1            TO  (-INFINITY,BOUND),
  36. C                     INF = 2             TO (-INFINITY,+INFINITY).
  37. C
  38. C            EPSABS - DOUBLE PRECISION
  39. C                     ABSOLUTE ACCURACY REQUESTED
  40. C            EPSREL - DOUBLE PRECISION
  41. C                     RELATIVE ACCURACY REQUESTED
  42. C                     IF  EPSABS.LE.0
  43. C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  44. C                     THE ROUTINE WILL END WITH IER = 6.
  45. C
  46. C            LIMIT  - INTEGER
  47. C                     GIVES AN UPPER BOUND ON THE NUMBER OF SUBINTERVALS
  48. C                     IN THE PARTITION OF (A,B), LIMIT.GE.1
  49. C
  50. C         ON RETURN
  51. C            RESULT - DOUBLE PRECISION
  52. C                     APPROXIMATION TO THE INTEGRAL
  53. C
  54. C            ABSERR - DOUBLE PRECISION
  55. C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
  56. C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
  57. C
  58. C            NEVAL  - INTEGER
  59. C                     NUMBER OF INTEGRAND EVALUATIONS
  60. C
  61. C            IER    - INTEGER
  62. C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
  63. C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
  64. C                             ACCURACY HAS BEEN ACHIEVED.
  65. C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE. THE
  66. C                             ESTIMATES FOR RESULT AND ERROR ARE LESS
  67. C                             RELIABLE. IT IS ASSUMED THAT THE REQUESTED
  68. C                             ACCURACY HAS NOT BEEN ACHIEVED.
  69. C                     IER.LT.0 EXIT REQUESTED FROM USER-SUPPLIED
  70. C                             FUNCTION.
  71. C
  72. C            ERROR MESSAGES
  73. C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
  74. C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE
  75. C                             SUBDIVISIONS BY INCREASING THE VALUE OF
  76. C                             LIMIT (AND TAKING THE ACCORDING DIMENSION
  77. C                             ADJUSTMENTS INTO ACCOUNT). HOWEVER,IF
  78. C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED
  79. C                             TO ANALYZE THE INTEGRAND IN ORDER TO
  80. C                             DETERMINE THE INTEGRATION DIFFICULTIES.
  81. C                             IF THE POSITION OF A LOCAL DIFFICULTY CAN
  82. C                             BE DETERMINED (E.G. SINGULARITY,
  83. C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
  84. C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
  85. C                             INTERVAL AT THIS POINT AND CALLING THE
  86. C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
  87. C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
  88. C                             SHOULD BE USED, WHICH IS DESIGNED FOR
  89. C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
  90. C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS
  91. C                             DETECTED, WHICH PREVENTS THE REQUESTED
  92. C                             TOLERANCE FROM BEING ACHIEVED.
  93. C                             THE ERROR MAY BE UNDER-ESTIMATED.
  94. C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS
  95. C                             AT SOME POINTS OF THE INTEGRATION
  96. C                             INTERVAL.
  97. C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
  98. C                             ROUNDOFF ERROR IS DETECTED IN THE
  99. C                             EXTRAPOLATION TABLE.
  100. C                             IT IS ASSUMED THAT THE REQUESTED TOLERANCE
  101. C                             CANNOT BE ACHIEVED, AND THAT THE RETURNED
  102. C                             RESULT IS THE BEST WHICH CAN BE OBTAINED.
  103. C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
  104. C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT
  105. C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
  106. C                             OF IER.
  107. C                         = 6 THE INPUT IS INVALID, BECAUSE
  108. C                             (EPSABS.LE.0 AND
  109. C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  110. C                             RESULT, ABSERR, NEVAL, LAST, RLIST(1),
  111. C                             ELIST(1) AND IORD(1) ARE SET TO ZERO.
  112. C                             ALIST(1) AND BLIST(1) ARE SET TO 0
  113. C                             AND 1 RESPECTIVELY.
  114. C
  115. C            ALIST  - DOUBLE PRECISION
  116. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  117. C                      LAST  ELEMENTS OF WHICH ARE THE LEFT
  118. C                     END POINTS OF THE SUBINTERVALS IN THE PARTITION
  119. C                     OF THE TRANSFORMED INTEGRATION RANGE (0,1).
  120. C
  121. C            BLIST  - DOUBLE PRECISION
  122. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  123. C                      LAST  ELEMENTS OF WHICH ARE THE RIGHT
  124. C                     END POINTS OF THE SUBINTERVALS IN THE PARTITION
  125. C                     OF THE TRANSFORMED INTEGRATION RANGE (0,1).
  126. C
  127. C            RLIST  - DOUBLE PRECISION
  128. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  129. C                      LAST  ELEMENTS OF WHICH ARE THE INTEGRAL
  130. C                     APPROXIMATIONS ON THE SUBINTERVALS
  131. C
  132. C            ELIST  - DOUBLE PRECISION
  133. C                     VECTOR OF DIMENSION AT LEAST LIMIT,  THE FIRST
  134. C                     LAST ELEMENTS OF WHICH ARE THE MODULI OF THE
  135. C                     ABSOLUTE ERROR ESTIMATES ON THE SUBINTERVALS
  136. C
  137. C            IORD   - INTEGER
  138. C                     VECTOR OF DIMENSION LIMIT, THE FIRST K
  139. C                     ELEMENTS OF WHICH ARE POINTERS TO THE
  140. C                     ERROR ESTIMATES OVER THE SUBINTERVALS,
  141. C                     SUCH THAT ELIST(IORD(1)), ..., ELIST(IORD(K))
  142. C                     FORM A DECREASING SEQUENCE, WITH K = LAST
  143. C                     IF LAST.LE.(LIMIT/2+2), AND K = LIMIT+1-LAST
  144. C                     OTHERWISE
  145. C
  146. C            LAST   - INTEGER
  147. C                     NUMBER OF SUBINTERVALS ACTUALLY PRODUCED
  148. C                     IN THE SUBDIVISION PROCESS
  149. C
  150. C***REFERENCES  (NONE)
  151. C***ROUTINES CALLED  D1MACH,DQELG,DQK15I,DQPSRT
  152. C***END PROLOGUE  DQAGIE
  153.       DOUBLE PRECISION ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,
  154.      *  A2,BLIST,BOUN,BOUND,B1,B2,CORREC,DABS,DEFABS,DEFAB1,DEFAB2,
  155.      *  DMAX1,DRES,D1MACH,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,ERLAST,
  156.      *  ERRBND,ERRMAX,ERROR1,ERROR2,ERRO12,ERRSUM,ERTEST,F,OFLOW,RESABS,
  157.      *  RESEPS,RESULT,RES3LA,RLIST,RLIST2,SMALL,UFLOW
  158.       INTEGER ID,IER,IERRO,INF,IORD,IROFF1,IROFF2,IROFF3,JUPBND,K,KSGN,
  159.      *  KTMIN,LAST,LIMIT,MAXERR,NEVAL,NRES,NRMAX,NUMRL2
  160.       LOGICAL EXTRAP,NOEXT
  161. C
  162.       DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT),
  163.      *  RES3LA(3),RLIST(LIMIT),RLIST2(52)
  164. C
  165.       EXTERNAL F
  166. C
  167. C            THE DIMENSION OF RLIST2 IS DETERMINED BY THE VALUE OF
  168. C            LIMEXP IN SUBROUTINE DQELG.
  169. C
  170. C
  171. C            LIST OF MAJOR VARIABLES
  172. C            -----------------------
  173. C
  174. C           ALIST     - LIST OF LEFT END POINTS OF ALL SUBINTERVALS
  175. C                       CONSIDERED UP TO NOW
  176. C           BLIST     - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS
  177. C                       CONSIDERED UP TO NOW
  178. C           RLIST(I)  - APPROXIMATION TO THE INTEGRAL OVER
  179. C                       (ALIST(I),BLIST(I))
  180. C           RLIST2    - ARRAY OF DIMENSION AT LEAST (LIMEXP+2),
  181. C                       CONTAINING THE PART OF THE EPSILON TABLE
  182. C                       WICH IS STILL NEEDED FOR FURTHER COMPUTATIONS
  183. C           ELIST(I)  - ERROR ESTIMATE APPLYING TO RLIST(I)
  184. C           MAXERR    - POINTER TO THE INTERVAL WITH LARGEST ERROR
  185. C                       ESTIMATE
  186. C           ERRMAX    - ELIST(MAXERR)
  187. C           ERLAST    - ERROR ON THE INTERVAL CURRENTLY SUBDIVIDED
  188. C                       (BEFORE THAT SUBDIVISION HAS TAKEN PLACE)
  189. C           AREA      - SUM OF THE INTEGRALS OVER THE SUBINTERVALS
  190. C           ERRSUM    - SUM OF THE ERRORS OVER THE SUBINTERVALS
  191. C           ERRBND    - REQUESTED ACCURACY MAX(EPSABS,EPSREL*
  192. C                       ABS(RESULT))
  193. C           *****1    - VARIABLE FOR THE LEFT SUBINTERVAL
  194. C           *****2    - VARIABLE FOR THE RIGHT SUBINTERVAL
  195. C           LAST      - INDEX FOR SUBDIVISION
  196. C           NRES      - NUMBER OF CALLS TO THE EXTRAPOLATION ROUTINE
  197. C           NUMRL2    - NUMBER OF ELEMENTS CURRENTLY IN RLIST2. IF AN
  198. C                       APPROPRIATE APPROXIMATION TO THE COMPOUNDED
  199. C                       INTEGRAL HAS BEEN OBTAINED, IT IS PUT IN
  200. C                       RLIST2(NUMRL2) AFTER NUMRL2 HAS BEEN INCREASED
  201. C                       BY ONE.
  202. C           SMALL     - LENGTH OF THE SMALLEST INTERVAL CONSIDERED UP
  203. C                       TO NOW, MULTIPLIED BY 1.5
  204. C           ERLARG    - SUM OF THE ERRORS OVER THE INTERVALS LARGER
  205. C                       THAN THE SMALLEST INTERVAL CONSIDERED UP TO NOW
  206. C           EXTRAP    - LOGICAL VARIABLE DENOTING THAT THE ROUTINE
  207. C                       IS ATTEMPTING TO PERFORM EXTRAPOLATION. I.E.
  208. C                       BEFORE SUBDIVIDING THE SMALLEST INTERVAL WE
  209. C                       TRY TO DECREASE THE VALUE OF ERLARG.
  210. C           NOEXT     - LOGICAL VARIABLE DENOTING THAT EXTRAPOLATION
  211. C                       IS NO LONGER ALLOWED (TRUE-VALUE)
  212. C
  213. C            MACHINE DEPENDENT CONSTANTS
  214. C            ---------------------------
  215. C
  216. C           EPMACH IS THE LARGEST RELATIVE SPACING.
  217. C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
  218. C           OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
  219. C
  220. C***FIRST EXECUTABLE STATEMENT  DQAGIE
  221.        EPMACH = D1MACH(4)
  222. C
  223. C           TEST ON VALIDITY OF PARAMETERS
  224. C           -----------------------------
  225. C
  226.       IER = 0
  227.       NEVAL = 0
  228.       LAST = 0
  229.       RESULT = 0.0D+00
  230.       ABSERR = 0.0D+00
  231.       ALIST(1) = 0.0D+00
  232.       BLIST(1) = 0.1D+01
  233.       RLIST(1) = 0.0D+00
  234.       ELIST(1) = 0.0D+00
  235.       IORD(1) = 0
  236.       IF(EPSABS.LE.0.0D+00.AND.EPSREL.LT.DMAX1(0.5D+02*EPMACH,0.5D-28))
  237.      *  IER = 6
  238.        IF(IER.EQ.6) GO TO 999
  239. C
  240. C
  241. C           FIRST APPROXIMATION TO THE INTEGRAL
  242. C           -----------------------------------
  243. C
  244. C           DETERMINE THE INTERVAL TO BE MAPPED ONTO (0,1).
  245. C           IF INF = 2 THE INTEGRAL IS COMPUTED AS I = I1+I2, WHERE
  246. C           I1 = INTEGRAL OF F OVER (-INFINITY,0),
  247. C           I2 = INTEGRAL OF F OVER (0,+INFINITY).
  248. C
  249.       BOUN = BOUND
  250.       IF(INF.EQ.2) BOUN = 0.0D+00
  251.       CALL DQK15I(F,BOUN,INF,0.0D+00,0.1D+01,RESULT,ABSERR,
  252.      *  DEFABS,RESABS,IER)
  253.       IF (IER .LT. 0) RETURN
  254. C
  255. C           TEST ON ACCURACY
  256. C
  257.       LAST = 1
  258.       RLIST(1) = RESULT
  259.       ELIST(1) = ABSERR
  260.       IORD(1) = 1
  261.       DRES = DABS(RESULT)
  262.       ERRBND = DMAX1(EPSABS,EPSREL*DRES)
  263.       IF(ABSERR.LE.1.0D+02*EPMACH*DEFABS.AND.ABSERR.GT.ERRBND) IER = 2
  264.       IF(LIMIT.EQ.1) IER = 1
  265.       IF(IER.NE.0.OR.(ABSERR.LE.ERRBND.AND.ABSERR.NE.RESABS).OR.
  266.      *  ABSERR.EQ.0.0D+00) GO TO 130
  267. C
  268. C           INITIALIZATION
  269. C           --------------
  270. C
  271.       UFLOW = D1MACH(1)
  272.       OFLOW = D1MACH(2)
  273.       RLIST2(1) = RESULT
  274.       ERRMAX = ABSERR
  275.       MAXERR = 1
  276.       AREA = RESULT
  277.       ERRSUM = ABSERR
  278.       ABSERR = OFLOW
  279.       NRMAX = 1
  280.       NRES = 0
  281.       KTMIN = 0
  282.       NUMRL2 = 2
  283.       EXTRAP = .FALSE.
  284.       NOEXT = .FALSE.
  285.       IERRO = 0
  286.       IROFF1 = 0
  287.       IROFF2 = 0
  288.       IROFF3 = 0
  289.       KSGN = -1
  290.       IF(DRES.GE.(0.1D+01-0.5D+02*EPMACH)*DEFABS) KSGN = 1
  291. C
  292. C           MAIN DO-LOOP
  293. C           ------------
  294. C
  295.       DO 90 LAST = 2,LIMIT
  296. C
  297. C           BISECT THE SUBINTERVAL WITH NRMAX-TH LARGEST ERROR ESTIMATE.
  298. C
  299.         A1 = ALIST(MAXERR)
  300.         B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR))
  301.         A2 = B1
  302.         B2 = BLIST(MAXERR)
  303.         ERLAST = ERRMAX
  304.         CALL DQK15I(F,BOUN,INF,A1,B1,AREA1,ERROR1,RESABS,DEFAB1,IER)
  305.         IF (IER .LT. 0) RETURN
  306.         CALL DQK15I(F,BOUN,INF,A2,B2,AREA2,ERROR2,RESABS,DEFAB2,IER)
  307.         IF (IER .LT. 0) RETURN
  308. C
  309. C           IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
  310. C           AND ERROR AND TEST FOR ACCURACY.
  311. C
  312.         AREA12 = AREA1+AREA2
  313.         ERRO12 = ERROR1+ERROR2
  314.         ERRSUM = ERRSUM+ERRO12-ERRMAX
  315.         AREA = AREA+AREA12-RLIST(MAXERR)
  316.         IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2)GO TO 15
  317.         IF(DABS(RLIST(MAXERR)-AREA12).GT.0.1D-04*DABS(AREA12)
  318.      *  .OR.ERRO12.LT.0.99D+00*ERRMAX) GO TO 10
  319.         IF(EXTRAP) IROFF2 = IROFF2+1
  320.         IF(.NOT.EXTRAP) IROFF1 = IROFF1+1
  321.    10   IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF3 = IROFF3+1
  322.    15   RLIST(MAXERR) = AREA1
  323.         RLIST(LAST) = AREA2
  324.         ERRBND = DMAX1(EPSABS,EPSREL*DABS(AREA))
  325. C
  326. C           TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
  327. C
  328.         IF(IROFF1+IROFF2.GE.10.OR.IROFF3.GE.20) IER = 2
  329.         IF(IROFF2.GE.5) IERRO = 3
  330. C
  331. C           SET ERROR FLAG IN THE CASE THAT THE NUMBER OF
  332. C           SUBINTERVALS EQUALS LIMIT.
  333. C
  334.         IF(LAST.EQ.LIMIT) IER = 1
  335. C
  336. C           SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
  337. C           AT SOME POINTS OF THE INTEGRATION RANGE.
  338. C
  339.         IF(DMAX1(DABS(A1),DABS(B2)).LE.(0.1D+01+0.1D+03*EPMACH)*
  340.      *  (DABS(A2)+0.1D+04*UFLOW)) IER = 4
  341. C
  342. C           APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
  343. C
  344.         IF(ERROR2.GT.ERROR1) GO TO 20
  345.         ALIST(LAST) = A2
  346.         BLIST(MAXERR) = B1
  347.         BLIST(LAST) = B2
  348.         ELIST(MAXERR) = ERROR1
  349.         ELIST(LAST) = ERROR2
  350.         GO TO 30
  351.    20   ALIST(MAXERR) = A2
  352.         ALIST(LAST) = A1
  353.         BLIST(LAST) = B1
  354.         RLIST(MAXERR) = AREA2
  355.         RLIST(LAST) = AREA1
  356.         ELIST(MAXERR) = ERROR2
  357.         ELIST(LAST) = ERROR1
  358. C
  359. C           CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING
  360. C           IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
  361. C           WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
  362. C
  363.    30   CALL DQPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX)
  364.         IF(ERRSUM.LE.ERRBND) GO TO 115
  365.         IF(IER.NE.0) GO TO 100
  366.         IF(LAST.EQ.2) GO TO 80
  367.         IF(NOEXT) GO TO 90
  368.         ERLARG = ERLARG-ERLAST
  369.         IF(DABS(B1-A1).GT.SMALL) ERLARG = ERLARG+ERRO12
  370.         IF(EXTRAP) GO TO 40
  371. C
  372. C           TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE
  373. C           SMALLEST INTERVAL.
  374. C
  375.         IF(DABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90
  376.         EXTRAP = .TRUE.
  377.         NRMAX = 2
  378.    40   IF(IERRO.EQ.3.OR.ERLARG.LE.ERTEST) GO TO 60
  379. C
  380. C           THE SMALLEST INTERVAL HAS THE LARGEST ERROR.
  381. C           BEFORE BISECTING DECREASE THE SUM OF THE ERRORS OVER THE
  382. C           LARGER INTERVALS (ERLARG) AND PERFORM EXTRAPOLATION.
  383. C
  384.         ID = NRMAX
  385.         JUPBND = LAST
  386.         IF(LAST.GT.(2+LIMIT/2)) JUPBND = LIMIT+3-LAST
  387.         DO 50 K = ID,JUPBND
  388.           MAXERR = IORD(NRMAX)
  389.           ERRMAX = ELIST(MAXERR)
  390.           IF(DABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90
  391.           NRMAX = NRMAX+1
  392.    50   CONTINUE
  393. C
  394. C           PERFORM EXTRAPOLATION.
  395. C
  396.    60   NUMRL2 = NUMRL2+1
  397.         RLIST2(NUMRL2) = AREA
  398.         CALL DQELG(NUMRL2,RLIST2,RESEPS,ABSEPS,RES3LA,NRES)
  399.         KTMIN = KTMIN+1
  400.         IF(KTMIN.GT.5.AND.ABSERR.LT.0.1D-02*ERRSUM) IER = 5
  401.         IF(ABSEPS.GE.ABSERR) GO TO 70
  402.         KTMIN = 0
  403.         ABSERR = ABSEPS
  404.         RESULT = RESEPS
  405.         CORREC = ERLARG
  406.         ERTEST = DMAX1(EPSABS,EPSREL*DABS(RESEPS))
  407.         IF(ABSERR.LE.ERTEST) GO TO 100
  408. C
  409. C            PREPARE BISECTION OF THE SMALLEST INTERVAL.
  410. C
  411.    70   IF(NUMRL2.EQ.1) NOEXT = .TRUE.
  412.         IF(IER.EQ.5) GO TO 100
  413.         MAXERR = IORD(1)
  414.         ERRMAX = ELIST(MAXERR)
  415.         NRMAX = 1
  416.         EXTRAP = .FALSE.
  417.         SMALL = SMALL*0.5D+00
  418.         ERLARG = ERRSUM
  419.         GO TO 90
  420.    80   SMALL = 0.375D+00
  421.         ERLARG = ERRSUM
  422.         ERTEST = ERRBND
  423.         RLIST2(2) = AREA
  424.    90 CONTINUE
  425. C
  426. C           SET FINAL RESULT AND ERROR ESTIMATE.
  427. C           ------------------------------------
  428. C
  429.   100 IF(ABSERR.EQ.OFLOW) GO TO 115
  430.       IF((IER+IERRO).EQ.0) GO TO 110
  431.       IF(IERRO.EQ.3) ABSERR = ABSERR+CORREC
  432.       IF(IER.EQ.0) IER = 3
  433.       IF(RESULT.NE.0.0D+00.AND.AREA.NE.0.0D+00)GO TO 105
  434.       IF(ABSERR.GT.ERRSUM)GO TO 115
  435.       IF(AREA.EQ.0.0D+00) GO TO 130
  436.       GO TO 110
  437.   105 IF(ABSERR/DABS(RESULT).GT.ERRSUM/DABS(AREA))GO TO 115
  438. C
  439. C           TEST ON DIVERGENCE
  440. C
  441.   110 IF(KSGN.EQ.(-1).AND.DMAX1(DABS(RESULT),DABS(AREA)).LE.
  442.      * DEFABS*0.1D-01) GO TO 130
  443.       IF(0.1D-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1D+03.
  444.      *OR.ERRSUM.GT.DABS(AREA)) IER = 6
  445.       GO TO 130
  446. C
  447. C           COMPUTE GLOBAL INTEGRAL SUM.
  448. C
  449.   115 RESULT = 0.0D+00
  450.       DO 120 K = 1,LAST
  451.         RESULT = RESULT+RLIST(K)
  452.   120 CONTINUE
  453.       ABSERR = ERRSUM
  454.   130 NEVAL = 30*LAST-15
  455.       IF(INF.EQ.2) NEVAL = 2*NEVAL
  456.       IF(IER.GT.2) IER=IER-1
  457.   999 RETURN
  458.       END
  459.